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Abstract 

In this paper, a lattice Boltzmann equation (LBE) model is proposed for binary fluids based on 
a quasi-incompressible phase-field model [J. Shen et al. Comm. Comp. Phys. 13, 1045 (2013)]. 
Compared with the other incompressible LBE models based on the incompressible phase-field 
theory, the quasi-incompressible model conserves mass locally. A series of numerical simulations 
are performed to validate the proposed model, and comparisons with an incompressible LBE model 
[H. Liang et al, Phys. Rev. E 89, 053320 (2014)] are also carried out. It is shown that the proposed 
model can track the interface accurately, and the predictions by the quasi-incompressible and 
incompressible models agree qualitatively well as the distribution of chemical potential is uniform, 
otherwise differ significantly. 
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I. INTRODUCTION 


Multiphase fluid flows are ubiquitous in engineering problems and natnral processes. 


Generally, a phase interface can be described by sharp interface approach [ll-y] or diffuse 
interface approach j4-ll|. In the sharp interface approach, the fluid is separated into some 
sub-domains by sharp interfaces in which each sub-domain contains only one phase, and the 
fluid properties such as density and viscosity are discontinnous across the interfaces. On 
the contrary, in diffuse interface approach, the fluid is treated continuously in the whole 
domain and the flnid properties vary smoothly across interfaces. An attractive feature 
of the diffuse interface method is it can model the complex interfacial dynamics without 
explicitly tracking the interfaces, and this feature makes it an ideal basis for developing 
efficient numerical schemes. 

In the diffuse interface approach, usually a phase-field variable (or order parameter) is 
used to distinguish different phases. The variable takes two distinct constant values in the 
bulk regions of the two phases, respectively, and changes smoothly across the interface. 
Based on the phase-field variable and its gradient, the free-energy of the system can be 
modelled, from which one can obtain a transport equation for the order parameter. The dy¬ 
namic change of the phase interface can then be described by this equation coupled with the 
governing eqnations of the flow. In most of previous works on immiscible binary mixtures 
of incompressible fluids, the flow is usually assumed to be governed by the incompressible 
Navier-Stokes equations including the interfacial force. However, as pointed out in j^, the 
assumption that the mixture is incompressible in the whole region is inconsistent with the 
conservation of mass as the densities of the fluids are unequal. To remedy this physical prob¬ 
lem, a quasi-incompressible phase-field model, which assumes the mixture is incompressible 


in bnlk regions bnt compressible in the mixing layer, has been proposed 


0 . 


A number of numerical schemes have been developed based on phase-field models inclnd- 






ing spectral methods (^, ll3, llSj , finite element methods |l4l-ll6l| and LBE methods |9l-lll| . 
Among these methods, the LBE method has received particular attentions due to some dis¬ 


tinctive featnres 


i. B 


which 


j^. The first phase-field LBE model was proposed by He et al. 
adopts an order parameter to track the interface of two incompressible fluids. However, there 
exist some differences between the derived governing equations and the phase-field theory 
for incompressible two-phase flows 23, and numerical instability can be produced for 
systems with a large density ratio. Later some improved LBE models based on phase-held 
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theory have been developed from different viewpoints. For instance, in order to improve 
numerical stability, Lee and Lin jn] designed a three-stage discretization multiphase lattice 
Boltzmann (LB) scheme by discretizing the gradient terms in different manners before and 
after the stea^iag step. Later, Fakhari et at. further ge^ralized thyucde. fl by 


employing a multi-relaxation-time collision operator. Zheng m and Zu (23[ respectively 
proposed the modified LBE models in order to recover the correct Cahn-Hilliard (CH) equa¬ 
tion. However, the extra terms in the both models will produce a large error in the interface 
capturing and the computation will become unstable as the dimensionless relaxation time 
equals to ij^. To overcome these problems, Liang et al. 24] proposed a new LBE model 
by introducing a time-dependent source term in the evolution equation. Recently, Zheng et 


al. 


251 ] presented an alternative model based on the kinetic theory to solve the problem. 


Although those LBE models 


0,y, 


25| could recover the CH equation exactly, the recov¬ 


ered momentum equations are still inconsistent with the target momentum equations for the 
incompressible flows. Li et al. [2^ noted this problem and proposed a correction method 
by introducing an artificial interfacial force. 

All of the above LBE models were developed based on the incompressible phase-field 
models that do not conserve mass locally as the two fluids have different densities. In this 
work we aim to develop a LBE model based on the quasi-incompressible phase-field theory 
for two-phase flows, which can ensure the exact mass conservation. The rest of this paper 
is organized as follows. In Sec. 2, the quasi-incompressible phase held model is briefly 
reviewed, and a LBE model based on this theory is constructed in Sec. 3. In Sec. 4, some 
numerical simulations are carried out to validate the proposed model, and some comparisons 
with a recent incompressible LBE model are also made. A brief summary is presented in 
section 5. 


II. QUASI-INCOMPRESSIBLE PHASE-FIELD MODEL 

In the phase held theory for a two-phase system, the thermodynamic behavior can be 
described by a free energy function with respect to an order parameter (p 

F(^)=/'[</>w + |lv0nda (1) 

J 

where p is used to distinguish different phases, is the bulk free-energy density, k is the 
coefficient of surface tension, and Q is the control volume. 
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For binary fluids, a double-well form of free-energy density 


0,El 


can be used 


^(0) = (^{4> - 


( 2 ) 


where and are the equilibrium values of the order parameters for fluids A and B, 


respectively, B is a constant related to the interfacial thickness W js, l28|, l29| and the surface 


tension a j29l . 


and 


m, 


w = 


ISk 


“ I V /3 ’ 


a = 




6 


( 3 ) 


( 4 ) 


With the bulk free energy, the chemical potential fi 

SF d4> 


B 


28|, 


29| can be obtained 


= 4 / 3(0 - 0 u )(0 - 0 b )(0 - 


( 5 ) 


and the order-parameter profile across the equilibrium interface can be obtained by solving 

K(f>)= 0 y , 

rh A — U 7-. A — /d^ T-. A 0/* \ 

( 6 ) 


( ) ^ ^ f 2C 

2 2 \W 


where 0 is the coordinate normal to the interface. T 


ion js 


be described by the Cahn-Hilliard (CH) equation 


d(f) 


le evolution of the order parameter can 


m, 


31 


^ + V-(0u) = V-(AVAi), 


( 7 ) 


where A is the mobility coefficient and u is the fluid velocity. 

In the incompressible phase-field model, the fluid is assumed to be incompressible every¬ 
where, and the flow can be described by the incompressible Navier-Stokes equations with 


an interfacial force 


B 


m, 


^ + V-(0u) = V-(AVAi), 


p -I- u • Vu^ = -Vp -I- V • (Vu Vu^)] F, 

V-u = 0, 


with 


P 


4>A — (t>B 


(Pa — Pb) + Pb, 


( 8 ) 

( 9 ) 

( 10 ) 

( 11 ) 
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where pa and ps are the densities of fluids A and B, respectively. From Eqs. (jH]), (fTIH) and 
dm, we can obtain 

f)n rl n 

( 12 ) 


| + V.(P») = |V.|AV.]. 


where 


dp Pa — Pb 


(13) 


dtp 4 >a — pB 

It is obvious that the mass conservation is constrained by the dp/dp and V ■ [AV/r]. In 
generally, V • [AV/u] is nonzero in the interfacial region. Hence, as long as pa ^ Pb, the 
mass is not locally conserved in the incompressible phase-field model. 


In the quasi-incompressible phase held model 


271], the governing equations are expressed 


as 


dp 

dt 


+ V • ((/)u) = V • (AV/i), 


P 


(?U 


(14) 

(15) 

(16) 

— . ( 17 ) 

p - <pdp/d4j 4>A - 4>BPr ’ 
where p is the hydrodynamic pressure, b is the kinematic viscosity, F is the total force 

including the surface tension force Fs(= —p'Vp) and other body forces F^, Pr is the density 
ratio Pa/Pb- From Eqs. flT^ . (ITHl) and (IMl) . one can obtain that 


with 


+ u ■ VuJ = —Vp + V ■ [piA (Vu + Vu^)] + F, 
V • u = —yV • [AV/i], 

dp/dp 


7 


Pr 


flll> + V ■ (/3U) = 0, 


(18) 


which means that the mass is conserved locally in the qusi-incompressible model. Further¬ 
more, equation (IT^ suggests that the huid is compressible in the mixing zone. To investigate 
the effect of compressibility, we substitute Eq. flTHl) into Eq. m to get, 

dp 


dt 


+ u ■ 'Vp = (1 -|- 70)V ■ (AV/r). 


(19) 


On the other hand, from Eqs. dHj) and 


dp 


we can obtain 


— + u-Vp = V ■ {XVp). 


( 20 ) 


From Eqs. flT^ and (EOl), we can see that the discrepancy between the two models is related 
to the term jpV ■ (AVp), which depends on the density ratio and the spatial distribution 
of the chemical potential. If the chemical potential is uniformly distributed, the additional 
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term plays weak role in the results, otherwise the discrepancy between the two models is 
tremendous. Eq. m can be also expressed in dimensionless formulation as 


dt(l> + u-V(l>=il + 7(/))Pe-iV2{4((/) - - (0n + 0s)/2] 

-C'n2(0^-0B)'V20/8}, 


( 21 ) 


where Pe = UcLc/XP is the Peclet number and Cn = W/Lc is the Cahn number with the 
characteristic length Lc and velocity Uc- This suggests that the distinctions between the 
two models are also related to the magnitudes of Pe and Cn. In the numerical simulations, 
we will investigate the difference between the two models by changing the dimensionless 
parameters Pe, Cn and 7 . 


III. THE QUASI-INCOMPRESSIBLE LBE MODEL 


In this section, we will propose the LBE model based on the quasi-incompressible phase- 


field equations 


27|. The model consists of two LBEs, one for the CH equation, and one for 


the Navier-Stokes equations, 


/*(x -h Ci5t, t + 5t)- t) = -[/i(x, t) - /^''(x, t)] + St 


1 - 


2Tf 


F„ ( 22 ) 


gi(x + dSt, t + 5t)- gi(x, t) = -[gi(x, t) - gf (x, t)] + St 

To- 


1 

2 tc 


gJ 


Gi, (23) 


where /j(x, t) and gj(x, t) are the distribution functions for the hydrodynamics and order 
parameter fields, respectively, Cj is the discrete velocity in the tth direction, St is the time 
step, Tf and Tg are dimensionless relaxation times related to the shear viscosity and mobility, 
respectively. Pi and Gi are the distribution functions for the force term, and Gi is used for 
eliminating the extra term in the CH equation jsj]. The local equilibrium distribution 
functions /®'^(x, t) and g^‘^{x,t) are respectively defined as 


/r= i^iip+4psi{u)], 


(24) 


gf = Hi + uJi(l)Si{u), 


(25) 


with 


sdu) = 


u uu : (cjCj — C3I 




(26) 
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0 — (1 — uo)aii, i = 0 
oJian, i 0 


(27) 


H,= 


where coi is the weighting coefficient, D is the spatial dimension and Cg is the sound speed 
for an ideal fluid, and a is an adjustable parameter for a given mobility. 

In Eqs. (l 2 ^ and (ESI), the source terms F* and Gi are respectively given by 


Fi = (cj - u) • [a;jFri(u) + u)iSi{u)clVp] - uJiclfqV ■ (AV/r), 

Gi = —^ (cj - u) • (Vp - F) a;iri(u). 

cjp 

The macroscopic quantities, 0, u and p, are computed evaluated as 


i 



i 

The kinetic viscosity u and the mobility A are respectively given by 


(28) 

(29) 

(30) 

(31) 

(32) 


u — cl{Tf — 0.5)5t, A = c1{Tg — 0.5)Q;(5t. 


(33) 


Through the Chapman-Enskog analysis (see the Appendix for details), we can obtain the 
following macroscopic hydrodynamics equations 

= —Vp + V • [pu (Vu + Vu^)] + F, (34) 

-L!^ + V ■ u = -tV ■ [AV/j], (35) 

cjpdt 

^ + V-(0u) = V-(AVp). (36) 

In the limit of low Mach number (Ma = |u|/cs), the dynamic pressure is assumed to be 
p ~ O(Ma^), and the above set of equations reduce to quasi-incompressible model given by 
Eqs. IHM to flT 6 ll . 

In the present work, we consider two-dimensional cases, and the two-dimensional nine- 
velocity (D2Q9) LBE model is used without loss of generality, in which Cq = (0, 0), Cj=i _4 = 
c{cos[(i — 1)7 t/2 ], sin[(i — l)7r/2]}, Cj= 5_8 = \/2c{cos[(2z — l)7r/4], sin[(2z — l)7r/4]}, and the 
corresponding weight coefficients are wq = 4/9, loi -4 = 1/9 and ws-s = 1/36. The sound 


P 


du 

'm 


+ u • Vu 
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speed Cg is given by Cg = where c = SxjSt, with 5x representing the lattice space. For 

simplicity, we set the lattice space and time increment as the length and time units, i.e., 
5x = 5t = 1. In the computations, the gradient operators are discretized with the isotropic 
central scheme y. 


IV. NUMERICAL SIMULATIONS 


In this section, we will validate the accuracy of the proposed quasi-incompressible LBM, 
and compare it with a recent LBE model based on the incompressible phase-field theory 


given m 


24| by a series of numerical simulations. 


A. One-dimensional flat interface 


We firstly validate the proposed LBE model by a flat interface test. Initially, the central 
region (25 < y < 75) is filled with fluid A and the rest is occupied by fluid B. The order 
parameter and density profiles are respectively set to be at equilibrium, i.e., 


+ ^^^tanhyi, y<50 
— ^a-<Pb tanh|/ 2 , y > 50 


(37) 


M^ + M^tanhi/i, 


Po{y) = < 


y<50 


(38) 


PA^ _ PA^ tanhi/ 2 , y>50 

where yi = 2{ij — 25) jW and 1/2 = 2(|/ — 75)/IF. The lattice used is x = 10 x 100 
and periodic boundary conditions are employed in both x and y directions. The other 
parameters are fixed as pa = Pb = 0.2, t/ = 1, = 1, = 1, = 0 and a = 0.001. 

The effects of Pe and Cn numbers on the density distributions are investigated. Note that 
the characteristic length and velocity in this paper take the values of lattice space 5x and 
velocity c, respectively. Figure [T] shows the density distribution across the interface for 
the quasi-incompressible (Quasi) model. From the figure, it can be seen that the density 
profiles match the analytical profile which indicates the accuracy of the proposed model in 
the interface tracking. 


8 














FIG. 1. (Color online) Density profiles across the interface with (a) different values of Pe with 
Cn = 4, and (b) different values of Cn with Pe = 1000. 


B. Stationary droplet 


A 2D stationary droplet problem is further tested to verify the present model. Initially, 
a circular droplet with radius ranging from 20 to 40 is placed in the middle of the compu¬ 
tational domain with Nx x Ny = 100 x 100. The initial order parameter and density fields 
profile are given by 


M^,y) = tank I 2^^ 


Po(x:, y) = o • 2 


sjix- Xcf + (?/ - ycY 



(39) 


(40) 


where (xcUc) is the center of the droplet. Different values of Pe and Cn are respectively 
investigated. The other parameters are fixed as pa = 1, Pb = 0.2, r/ = 1, Tg = 1, (f)A = 1, 
(pB — 0 and a = 0.001. When the droplet reaches the equilibrium state, the pressure 
difference AP between the inside and outside droplet should satisfy the Laplace law, i.e., 
AP = a/R, where P is calculated by P = po — + k|V0P/2 +p with the equation of 


state po = 4)d(j,'ip — ip (23|, l25|. Therefore, the surface tension can be calculated by (Jlbm = 


RAP, and the numerical predictions and theoretical values of the surface tension are shown 
in Table HI It can be seen that the present quasi-incompressible LBE model satisfies the 
Laplace law. The distributions of the order parameter predicted by the both LBE models 
are shown in Fig. El and no obvious difference can be observed. In order to observe the 
distinctions between them, a moving interface problem is attempted to be simulated in the 
next section. 
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TABLE 1. Numerical and theoretical values of surface tension with different Pe and Cn. 


Surface tension 

Pe {Cn = 4) 

Cn {Pe = 200) 

50 125 200 

4 6 8 

Numerical (xl0“^) 

0.992 0.992 0.992 

0.992 0.993 0.989 

Theoretical (xl0“^) 

1 1 1 

1 1 1 


Quasi-incompressible Incompressible 



Cn = 4 Cn = 4 


(/^ = 0.5 



20 40 60 80 100 


Cn = A 



Cn = 8 Cn = 8 


100 

80 

60 

40 

20 

20 40 60 80 100 

Cn = 8 



FIG. 2. (Color online) Order parameter configurations under different values of Cn with Pe = 200 
ffxed for the quasi-incompressible (Quasi) and incompressible (Incom) LBE models. 


C. Bubble rising under buoyancy 

In this section, a bubble rising under buoyancy is simulated to compare the two LBE 
models. Initially, a light circular bubble (fluid B) with radius R is immersed in another fluid 
(A) with higher density. To generate the buoyancy effect, a body force, Fh^y = —{p — pA)g, 
is added to the fluid flow, where g is the gravitational acceleration. In the simulations, the 
computational domain is set to be 160 x 480, and periodic boundary conditions are applied to 
all boundaries. The other parameters are set as follows: pa = 1, Pb = 0.5, g = 10“^, (f)A = 1, 
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4>b = 0, t/ = 1, Tg = 1, (t = 0.001, R = 32, Pe = 50 and Cn = 4. Figure [31 shows the shape 
of the rising bubble at different times predicted by the two LBE models. It can be seen that 
the results are quite similar. A comparison of the interface shapes at t = 10^ and 4 x 10^ 
confirms the similarity in Fig. jH Figure E] shows the distributions of the dynamic pressure 
at different times, and the difference in the vicinity of the interface is more obvious. Figures 
initoIHlshow the bubble velocity and the normalized velocity differences between the two LBE 
models. It can be seen from Eigs. (HlandlTl the horizontal and vertical velocity components 
are nearly identical, but from Eig. |Hl we can observe some the normalized velocity differences 
with maximum magnitude of order 10“^. Based on the above observations, we can conclude 
that the predictions by the two LBE models yield almost the same results for this test case. 
According to the previous theoretical analysis, these phenomena are reasonable since the 
initial equilibrium order parameter yields the approximate uniform chemical potential so 
that the term 70 V ■ (AV/r) exerts a weak influence on the results. 




(b)Incompressible 

FIG. 3. (Color online) Density confignration of the rising bubble at t/1000 = 10,20,30,35,40 for 
the quasi-incompressible model (a) and incompressible model (b). 
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FIG. 4. (Color online) Density profiles predicted by the quasi-incompressible (Quais) and incom¬ 
pressible (Incom) LBE models at t/1000 = 30,40. 




(a) Quasi-incompressible 



(b)Incompressible 

FIG. 5. (Color online) Dynamic pressure field of the rising bubble at t/1000 = 10, 20, 30, 35,40 for 
the quasi-incompressible model (a) and incompressible model (b). 
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(a) Quasi-incompressible 



(b)Incompressible 

FIG. 6. (Color online) Horizontal velocity of the rising bubble at t/1000 = 10, 20, 30, 35,40 for the 
quasi-incompressible model (a) and incompressible model (b). 
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(b)Incompressible 

FIG. 7. (Color online) Vertical velocity of the rising bubble at t/1000 = 10,20,30,35,40 for the 
quasi-incompressible model (a) and incompressible model (b). 
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(a) (b) 


FIG. 8. Normalized velocity difference between the quasi-incompressible (Quais) model and 
incompressible (Incom) model at t = 3 x 10^, i.e, (a) (itQuasi - ^tIncom)/^tmax,Quasi and (b) 
(^^Quasi “ ^^incom)/'ymax,Quasi) where Ui and Vj are the horizontal and vertical velocities of model 
i, respectively, and rtmax,Quasi and fmax,Quasi are the maximum horizontal and vertical velocities of 
quasi-incompressible model, respectively. 


D. Phase separation 

In this subsection, the simulation of phase separation will be carried out to further com¬ 
pare the two LBE models. Initially, the order parameter with a small perturbation is set as 
(p = 0.5[1 -h 0 . 1 sin( 47 ra:/La:)cos( 47 r?//Lr/)], where x, y represent the Cartesian coordinates, 
and Lx and Ly are the length and width of the computational domain, respectively. To 
consider the effects of Pe, Cn and 7 on the phase separation, different mobilities (A/3 = 0.02 
and 0.005), interfacial thicknesses {W = 4 and 8 ) and density ratios = 2 and 5) are 
taken into account in the simulations. The computational domain is set to be 100 x 100. 
The other parameters are set as pa = 1, 0a = 1) (pB = 0, r/ = 1, r^, = 1 and a = 0.001. 
Periodic boundary conditions are applied to the all boundaries. 

Figures [HI to [IS] depict the density, dynamic pressure and velocity of the mixture with 
various dimensionless parameters 7 , Pe and Cn. Firstly, we investigate the effect of 7 as 
shown in Figs. [HltodSl As 7 = 1 (see Figs. [Sjtodl]), in the transient period, the density, 
pressure and velocity fields predicted by the two models are quite different. In the steady 
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state, although the density helds appears to be similar, the dynamic pressure and velocities 
are remarkably distinct. As 7 = 4 (see Fig. US]), the density helds predicted by the two LBE 
models are opposite in the steady state, and the phase separation predicted by the present 
LBE model occurs earlier than that by the incompressible model. By comparing the Figs. 
[HI and [131 it can also be found that the density held from the present model varies with 7 
while that from the incompressible model does not. Then we further consider the ehect of 
Peclet number. As Pe increases to 200 (see Fig. m, it can be found the phase separation 
processes predicted by the two LBE models are both slowed down, but the distributions 
of the density helds are in opposite in the steady state. The effect of Cn on the phase 
separation process is also investigated. As shown in Fig. (131 when the Cn is increased from 
4 to 8, the density helds change greatly. In the transient period {t = 10^ to 9 x 10^), the 
results predicted by the two LBE models present similar conhgurations; while in the steady 
state, some strips with different angles of inclination appear in the density helds. The above 
phenomena completely exhibit the discrepancy between the two models when the chemical 
potential is nonuniformly distributed due to the non-equilibrium order parameter, and the 
discrepancy is deeply inhuenced by the dimensionless parameters 7, Pe and Cn. 




o°c 
o 



(a) Quasi-incompressible 


(b)Incompressible 


FIG. 9. (Color online) Density configuration with the minimum in blue and the maximum in red 
at t/1000 = 0.1,3,5,7,20 for the quasi-incompressible model (a) and incompressible model (b) as 
7 = 1 , Pe = 50, and Cn = 4. 
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(a) Quasi-incompressible 



(b)Incompressible 

FIG. 10. (Color online) Dynamic pressure field with the minimum in blue and the maximum in 
red at t/1000 = 0.1, 3, 5, 7, 20 for the quasi-incompressible model (a) and incompressible model (b) 
as 7 = 1, Pe = 50, and Cn = 4. 
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(b)Incompressible 


FIG. 11. (Color online) Horizontal velocity with the minimum in blue and the maximum in red 
at t/1000 = 0.1,3,5,7,20 for the quasi-incompressible model (a) and incompressible model (b) as 
7 = 1, Pe = 50, and Cn = 4. 
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(b)Incompressible 


FIG. 12. (Color online) Vertical velocity with the minimum in blue and the maximum in red at 
t/1000 = 0.1,3,5,7,20 for the quasi-incompressible model (a) and incompressible model (b) as 
7 = 1, Pe = 50, and Cn = 4. 
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(a) Quasi-incompressible 



(b)Incompressible 


FIG. 13. (Color online) Density configuration with the minimum in blue and the maximum in red 
at t/1000 = 0.1,1,2,4,15 for the quasi-incompressible model (a) and incompressible model (b) as 
7 = 4, Pe = 50, and Cn = 4. 
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(a) Quasi-incompressible 


(b)Incompressible 


FIG. 14. (Color online) Density configuration with the minimum in blue and the maximum in red 
at t/1000 = 0.1,10,20,25,50 for the quasi-incompressible model (a) and incompressible model (b) 
as 7 = 4, Pe = 200, and Cn = 4. 





(a) Quasi-incompressible 



(b)Incompressible 

FIG. 15. (Color online) Density configuration with the minimum in blue and the maximum in red 
at t/1000 = 10,60,90,200,240 for the quasi-incompressible model (a) and incompressible model 
(b) as 7 = 4, Pe = 200, and Cn — 8. 


V. CONCLUSIONS 

In this study, a LBE model for binary fluids is proposed based on the quasi-incompressible 
phase-field theory, which overcomes the mass conservation problem in the incompressible 
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phase-field LBM. To validate the accuracy of the proposed model and compare its perfor¬ 
mance with an incompressible model, a series of numerical tests are performed. 

Firstly, with the one-dimensional fiat interface and stationary droplet tests, it is shown 
that the proposed LBE model can track the interface accurately and satisfies the Laplace 
law. Furthermore, the comparison of the stationary droplet shows no obvious difference 
between the two LBE model. Then, the test of the bubble rising under buoyancy shows that 
some subtle distinctions exist in the two LBE models but overall they still agree with each 
other qualitatively. The results of the phase separation problem show the dynamic processes 
and static structures from the two models can be rather different. These phenomena indicate 
that there are some distinctions between the two models involved with the multiphase flows, 
especially for problems with the nonuniform distribution of chemical potential. Since the 
quasi-incompressible LBE model satisfies the fundamental mass conservation law, its results 
should be more reliable. 
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APPENDIX: CHAPMAN-ENSKOG ANALYSIS OF THE QUASI-INCOMPRESSIBLE 
LBE MODEL 

In this section, the proposed LBE model for hydrodynamic equations is firstly analyzed 
by applying the Chapman-Enskog expansion 

/i = ii‘°’ + T‘’ + £Vf+ ••■ («) 

d, = ed,,+e^d,„ V = £V„, Ft = ePf' + (42) 

with 

= - u) • ^F,F + uj,s,cl^p\ , fP= - n;,c2p7V.(AVAr), (43) 

where £ is a small expansion parameter. Using the Taylor expansion in Eq. (|22D . one can 
obtain 

A/i + f A"/. = -T/. - ft") + (' - ^) 
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where Tc = TfSt, and Di = dt + Ci ■ 'V. The substitution of Eqs. (IITI) and into Eq. (Ilil) 
yields the Chapman-Enskog system as 


0{e“) : f!°> = f!\ (45) 

0{e') : -Doi//"' - (l - ^) C = (46) 

O(e^): S..ii'"’+£>t>.C + |£'«/r = -^C+(4T) 

Then, the substitution of Eq. (Hbl) into (H7j) yields 

Meanwhile, from the definitions fjMl) and fH3l) . it is easy to calculate the following moments: 

i i 

XI CQ/r = ^Ip + c^puu, (50) 

i 

^ ^ ^/3^a7 “1“ 

i 

^i^f>=c;u.Vp, X^i5(‘) = _c>V.(AV,0. (52) 

X]c,Ff> = c=F, ^c.F.<‘'=0, (53) 

z z 

X CcFf = (F'u + uF') + c^u.Vp, X = - c^P7V.(AVp), (54) 

z z 

where F' = F + c,Wp. From Eqs. (El]), (E21) and (E^ . we can obtain 

E f'P = -fE C =f c;7pV.(AV,.). E C’ = 0, A,' > 2 (56) 

i i i 

= -^c^F, = 0. k> I (56) 

Z Z 

Taking the zeroth- and first-order moments of Eq. (EHl) , we can obtain 

-^dtoP -F V • u = 0, (57) 

pci 

dto (pu) -F V • (P + puu) = F. (58) 

Likewise, taking the zeroth- and first-order moments of Eq. fl48li . we can obtain 

d^P = -clfTfV • (AVp), (59) 
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dt, (c^pu) + V, 


0 • 


e(l) 


5t 


Vn- 


1 - 


2t, 




According to the Eqs. (IlHl) . and (IH7I) to 

- A ^ =clu ■ Vp + c^(uF' + F'u) + 4p(Vu + Vu”^) 

i 

-E(i- 2 ^) c.CjFf + 0(M»). 

Substituting Eq. flHTD into Eq. (EOI), we can obtain 

ai,(pu) = V - HVu + Vu^)], 


(0) 


= 0. (60) 


(61) 


(62) 


where u = Cg(r/ — 0.5)ht. From Eqs. (IHTIi and (IH^ . we can obtain the continuity equation 

1 


pcj 


dtp + V • u = —qV • (AVp). 


Similarly, the momentum equation can be derived from Eqs. fjSSll and fj621) 
dt{pn} + V ■ (puu) = —Vp + V • [pz/(Vu + Vu^)] + F. 


(63) 


(64) 


Next we will derive the CH equation from Eq. f|23l) by the Chapman-Enskog expansion. 
Similarly, the mnltiscale expansions are given by 


( 0 ) , ( 1 ) , 2 ( 2 ) , 

Si = Si ' +egy +8 gl + ..., 

dt = sdtQ + 8‘^dt^^ V = 6:Vo, Gi = 8G[\ 

Using the Taylor expansion in Eq. fj23l) . one can obtain 

Agj + = ““(Si “ SD A Gi, 


(65) 

( 66 ) 

(67) 


where Tc = r^St. Substituting Eqs. flHHl) and (IHHl) into the Eq. (IH7I) . we can obtain the 
following infinite consecutive series of equations 


0(£°) : gf^ = gf, 
0(£i): Doigf-fl- 


j-) = -i 

2tJ * r. 




0{e^) : dt^gf + Doigf^ + ^^gf^ = --gf^ 

2 Tc 


Then substituting Eq. into Eq. (1701) . we can obtain 


SiiSi + 1 “ y) + (y ’') = -4gf’. 


( 68 ) 

(69) 

(70) 

(71) 
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In order to recover the CH equation, we derive the following moments from the Eq. f|25l) and 

Eq. (I2nD, 

i i i 

^^Gj = 0, ^^CjGj = —G, (73) 

i i ^ 

where G = — Vp + F. Taking the zeroth-order moment of Eqs. (E!?D and (TTTI) . we can obtain 


9*00 + V ■ (0u) = 0, 


dtA 


St 


- Tr 


92 0 + 29*0 V • (</>u) + V • V • (0UU + clafi) - V • (^G) 

P . 


= 0 . 


According to Eqs. (1^ and fTTlD . the Eq. (TTSD reduces to the following equation 


(74) 

(75) 


9*10 = AV^/u, 


(76) 


where A = c^Q;5t(rg — 1/2). Combined Eqs. (1741) and flTHl) . we can obtain the CH equation 
as follows 

9*0 + V • (0u) = AVV- (77) 

From Eqs. flTTD . (1S41) and (1771) . we can derive the following mass conservation equation by 
neglecting the term 9*p, which is of order Ma^ 


9tp + V • (pu) = 0, 


(78) 


and then the momentum equation in Eq. (IMj) could be rewritten as 


9u 


p [ — + u • Vu I = —Vp + V • \^pp (Vu + Vu^)] + F. 


(79) 
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